Time propagation of constrained coupled Gaussian wave packets 
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The dynamics of quantum systems can be approximated by the time propagation of Gaussian wave 
packets. Applying a time dependent variational principle, the time evolution of the parameters of 
the coupled Gaussian wave packets can be calculated from a set of ordinary differential equations. 
Unfortunately, the set of equations is ill-behaved in most practical applications, depending on the 
number of propagated Gaussian wave packets, and methods for regularization are needed. We 
present a general method for regularization based on applying adequate nonholonomic inequality 
constraints to the evolution of the parameters, keeping the equations of motion well-behaved. The 
power of the method is demonstrated for a non-integrable system with two degrees of freedom. 

PACS numbers: 03.65.-w, 04.30.Nk 



I. INTRODUCTION 

The method of Gaussian wave packet propagation 
is a popular tool for quantum dynamics computations. 
Within this approximation it is assumed that an ini- 
tially Gaussian wave packet (GWP) stays Gaussian for 
all times. The time evolution of the wave packet is given 
by the time evolution of its parameters like width, phase, 
center, and momentum [ll]. For a single GWP, this rather 
crude approximation is in general only valid for short 
time propagation. The approximation can be signifi- 
cantly improved, if a superposition of GWP is used and 
these GWP are propagated in concert, since the num- 
ber of adjustable parameters is increased and the overall 
wave function is no longer restricted to a Gaussian shape 
0, S, 0, 0] ■ The equations of motion for the Gaussian pa- 
rameters are obtained from a time dependent variational 
principle (TDVP). It is well known that these coupled 
equations of motion for the time dependent parameters 
become ill-conditioned from time to time during the in- 
tegration depending on how many GWP are used. The 
reasons for the ill-conditioned behavior of the differential 
equations are near singularities of a matrix that has to be 
inverted after each time step of integration [1, |3) S S 13 ■ 
Using step size control the time steps of the integra- 
tion algorithm can become extremely small making the 
method impracticably slow. In the worst case even a 
failure of the numerical matrix inversion or the further 
integration may occur. 

Different solutions to this numerical problem were pro- 
posed, e.g. a regularization based on a singular value de- 
composition 0]. The singular value decomposition is ca- 
pable of regularizing the equations of motion in the sense 
that the method does not break down, however it does 
not solve the problem with the tiny step sizes Q. An- 
other proposal is to adjust the number of GWP during 
evolution by increasing or reducing their number depend- 
ing on whether the wave function spreads or shrinks to 
avoid redundancy [1, 0, 111 ■ 

It has also been discussed to simplify the equations of 
motion by keeping the widths of the propagated GWP 
fixed, called frozen Gaussian approximation 0, [^, [loj , 



or much cruder, to neglect the coupling between the 
GWP 3, 5]. Another proposal is to reduce the variational 
freedom by forcing the GWP to run on their classical tra- 
jectories 0, @, But of course these grave restrictions 
severely reduce the accuracy of the GWP method. 

Here we present a novel method to overcome the nu- 
merical problems or more precisely a method that avoids 
numerical problems in the first place. The idea is to im- 
pose adequate nonholonomic inequality constraints to the 
motion of each GWP, keeping the matrix regular. These 
constraints only become active when it is numerically 
necessary and otherwise leave the full variational freedom 
of the trial function. The method presented here is gen- 
eral and allows for the application of arbitrary (inequal- 
ity) constraints not only on GWP trial functions. There 
is numerical evidence, that near matrix singularities usu- 
ally result from widely varying amplitudes of largely over- 
lapping GWP. In our calculations it was sufficient to ac- 
count for one ingredient of the matrix singularity only, i.e. 
to constrain the amplitudes of the individual GWP to a 
reasonable domain. We account for the constraints in the 
time dependent variational principle and obtain different 
equations of motion as compared to the unconstrained 
variation. However, the equations of motion still have 
the form of a matrix equation as in the unconstrained 
case. Properly chosen constraints only slightly decrease 
the accuracy of the variational approximation. The addi- 
tional error introduced by the constraints decreases with 
a growing number of GWP. The method is able to avoid 
numerical problems rendering the integration by orders 
of magnitude faster. 

The article is organized as follows. In section HIl we re- 
capitulate the time dependent variational principle. The 
equations of motion for the Gaussian parameters ob- 
tained from the TDVP applied to GWP are given for 
completeness. In section Hill we account for the inequal- 
ity constraints in the TDVP and derive the regularized 
equations of motion. In section lTVl wc compare numerical 
results obtained from the GWP method with and without 
constraints in a two-dimensional non-integrable model 
potential, namely the 2D diamagnetic hydrogen atom. 
The accuracy of the constrained method is demonstrated 
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by comparison with other propagation techniques. A 
summary is given in section [V] 



II. TIME DEPENDENT VARIATIONAL 
PRINCIPLE 

The evolution of a quantum mechanical wave function 
is determined by the Schrodinger equation 

iiP{t) = Hi;{t) 

where the wave function V'(i) is an element of the Hilbert 
space. An approximate solution x(t) on a given manifold 
in Hilbert space can be obtained by a TDVP [13, fisl. [T3 . 
[isj . Here we choose the formulation of McLachlan [l4| . 
or equivalently the minimum error method Q , where the 
norm of the deviation between the right and the left hand 
side of the Schrodinger equation with respect to the trial 
function is to be minimized. The quantity 

J=||*0(O-i?xWll' = min 

is to be varied with respect to </) only, and then x = 4' 
is chosen. We assume the approximation manifold to 
be parametrized by a set of time dependent parameters 
z{t) = (zi(t),...,z„^(t)), i.e. xit) =x(zW)- In terms of 
these parameters the quantity / reads 
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which is a quadratic function of z for fixed values of z. 
The variation dcj) carries over to variations dz leading to 
the condition 



dl 



(2) 



',jr + izji one has the free- 
and dl/dzji = or to 



For complex parameters Zj = 
dom to take either dl/dzjr 
treat i* and Zj formally as independent parameters and 
to take either dl/dzj = or dl /dz* — 0. The resulting 
equations of motion are equivalent and read 
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in case of complex parameters z, where 
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The Hermitian matrix K is positive semi-definite since 
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> 0, 
(5) 



V c £ C"p , ensuring that the extremum of the quadratic 
quantity / is a minimum. 

The Schrodinger equation is replaced by a system of or- 
dinary first order differential equations of motion for the 
parameters z{t) where after every time step of integra- 
tion the set of simultaneous linear equations ([3]) must be 
solved for the time derivatives z if a numerical algorithm 
for ordinary differential equations, e.g. Runge-Kutta or 
Adams, is used. 



A. Application of the TDVP to GWP 

In this article a superposition of GWP as trial function 
is discussed. Each GWP (x e M-°) is of the form 



g(y'=,x) = ^i{{^-q!')A''{^-q!')+p''-{->^-q!')+i'' 



(6) 



where is a complex symmetric D x D matrix, the 
momenta p*^ and centers q*^ are real, D-dimensional vec- 
tors, and the phase and normalization are given by the 
complex scalars 7*^. The Gaussian parameters of the k-th 
GWP are denoted by y*^ — (A'^, p'^, q*^, 7*^). Their time 
argument is omitted for brevity. The trial function is a 
superposition of N such GWP 

N 

X(z,x)-^5(y^x), z=(yi,...,y^). (7) 

k=l 

Using a splitting of the Hamiltonian H = T+V we obtain 



^ 1 
Tx = 5]5(y',x)([*trA'=-7^ + p'=.(q^--p'=)] 



+ [-p^- + 2A'=(q^-p'=)]-(x-q^) 
-f(x-q''0[-i'=-2(A'=f](x-q'=)) 

^ 1 
^ E(^'o + • X + -xK,'=x) .g(y^ x), (8) 

k=l 

which defines, after sorting by powers of x, the complex 
scalars Vq, the complex vectors G and the com- 
plex symmetric D x D matrices V2 as the coefficients 
of a second order polynomial. According to the TDVP 
these coefficients (wq , vj;', Vj*^), k = 1,...,N are calcu- 
lated from a set of linear equations 



N 



N 



fe=l 

N 



k=l 



N N 

-E(5Vr^7x^2'x|/) =E(ffVr^"^(x)|/-); (9) 



fc=l 



fe=l 



; = l,...,iV; TO + n = 0,1,2; i,j = l,. 



On the right hand side the potential V{k) of the Hamil- 
tonian is inserted. It is straightforward to calculate the 
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time derivatives of the Gaussian parameters once the Un- 
ear equations (O are solved, since the differential equa- 
tions for the Gaussian parameters can be expressed by 
{vq, vJ, ^2*^), k = 1, . . . ,N according to their definition 
in equation ([8]): 

# = -2(yl'=)2 - iK^fe, 

p*^ = 2RcA'=s'= -Rcv^ -ReV^a'^q'': (10) 
^fc _ + itr A'^ + ^{p'')'^ — v'l ■ q'^ 

where s*^ — i(Im A'^)^^(Im vj;' + Im Vj'^q'^). Numerically 
it is more appropriate to introduce two additional D x D 
complex matrices B'', according to A'' = ^B'^{C'')^^ , 
and to integrate the equations of motion 

^ ~ ' fill 

instead of integrating A'' {t) directly, because the oscillat- 
ing {A''{t))'^ term causes numerical difficulties [lB|- For 
numerical accuracy, it is appropriate to symmetrize the 
matrix A''{t) after each time step. 

Equation ^ can be abbreviated by Kv — r when all 
coefficients {vQj-Vi^Vi), k = 1,...,N are put together 
into the complex vector v. All inner products in Hilbert 
space denoted by (.|.) are calculated in position space 
representation. The integrals that build up the compo- 
nents of the matrix K on the left hand side of equation 
© as well as the integrals on the right hand side can be 
solved analytically, provided the potential is of special 
form, e.g. polynomial, Gaussian or exponential. 

Given some initial wave function, i.e. the initial pa- 
rameters z{t = 0), the wave function is propagated by 
integrating the trajectories of the parameters. At every 
time step equation ([9]) must be solved for the coefficients 
V which are inserted in PH)) to obtain z. In the course 
of integration, depending on how many GWP are prop- 
agated in common, it will sooner or later happen that 
the matrix K associated with the set of linear equations 
becomes ill-conditioned, or even numerically singu- 
lar. As a result the time step of the integration routine 
becomes extremely small, rendering the method of GWP 
propagation impracticably slow. In the worst case, fur- 
ther integration or matrix inversion respectively, can even 
fail. 



III. INEQUALITY CONSTRAINED TDVP 

Matrix singularity problems arise from overcrowding 
the basis set, i.e. from situations where fewer GWP would 
be sufficient to represent the wave function. On the other 
hand for an accurate approximation of the wave function 
it is desirable to have a large number of adjustable pa- 
rameters. However, there is a discrepancy between the 



number of GWP necessary to give accurate results and 
the maximum number of GWP that can be propagated 
using the TDVP without numerical difficulties i. As 
mentioned above there exist different proposals to over- 
come this numerical problem, such as a singular value 
decomposition of the matrix Q or reducingthe num- 
ber of GWP when overcrowding takes place [1,13, Q. Also 
reducing the variational freedom by freezing the widths 
d, H, [a and choosing classical trajectories for the 
centers of the GWP [3, @,llH has been discussed. 

Our approach of regularizing the equations of motion 
for the parameters is based on minimizing the quantity 
/ in ll]) while certain inequality constraints are applied. 
The constraints must be chosen in such a way that they 
prevent the matrix iiT in ([9|) to become ill-conditioned. 
This means all Gaussian parameters evolve freely accord- 
ing to the TDVP, and the constraints only become active 
from time to time whenever the unconstrained evolution 
would drive the parameters in domains where the matrix 
would be too singular, and are switched off as soon as 
these 'forbidden' domains are left again. Formally spo- 
ken we reduce the space of admissible configurations to 
regions where the associated matrix K is regular. 

To demonstrate the generality of our method we first 
apply constraints to the general case of an arbitrary trial 
function whose parameters z{t) evolve according 

to equation ([3]). We derive their modified equations of 
motion which are obtained if the parameters z{t) are sub- 
ject to some arbitrary inequality constraints. Then we 
return to GWP trial functions (O and derive the mod- 
ification of equation ^ obtained when the GWP are 
subject to inequality constraints. Adequate constraints 
which prevent the matrix from singularity are presented 
and applied. 

Due to real inequality constraints it is convenient to use 
a real formulation of the equations. Complex quantities 
are split into their real and imaginary parts, which are 
denoted by the subscripts r and i, respectively. 

A. Inequality constrained TDVP on arbitrary trial 
functions 

Consider an arbitrary trial function and assume 

a real inequality constraint on the parameters z{t) G 0"^ 
which can be written in the form 

/(z,z*) = /(z,„zO = /(z) (12) 

where the function / is explicitly known. For brevity, the 
notation z = (z^, z^) G M^"j' will be used. 

As long as f{zr,Zi) > /min, all parameters evolve ac- 
cording to equation ^ without being affected by the re- 
striction. When f{zr, Zi) — /min is reached at some point 
in time t, the constraint becomes active, and we have to 
demand f{t) > 0, otherwise f{t + At) with some small 
positive At would violate the constraint ([H]). Therefore 
the quantity / of equation ([1]) at fixed z must be mini- 
mized with respect to z, where (zr,Zi) are now subject 
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to the constraint 

; df . ^df 

OZ,r OZi 



dl 
dz 



^ • z > 0. 



(13) 



In other words the possibly nonlinear constraint (|12p on z 
has been reduced to the linear constraint on z when 
/ = /min- Then the allowed domain of (z^, Zi) for search- 
ing the minimum of / is no more the whole space 
but the half-space / > linearly restricted by equation 
P3p . In general, minimization of a function on a given 
domain requires two steps, firstly to find the local inter- 
nal minima and secondly, to find the local minima on the 
boundaries. The global minimum in the given domain is 
obtained by comparison. Here it is sufficient to search for 
the minimum of / solely on the boundary of the domain 
defined by equation (llSp where the equality sign is ful- 
filled. That means the inequality (|13p may be replaced 
by the computationally much more feasible constraint 

— . z, + — • z, = — • z - 0. (14) 

The reason is that / is a positive definite parabolic func- 
tion of z whose absolute minimum lies outside the allowed 
domain by assumption. Since there are no internal min- 
ima / obviously takes its allowed minimum on the bound- 
ary of the allowed domain. The constraint is switched off 
again as soon as the trajectory z(t) of the absolute min- 
imum of / crosses the plane given by equation (jl4p in 
the (zr, Zi)-space at fixed values of (zj.,Zi). Note that 
arbitrary nonlinear constraints (jl2[) on z always lead to 
linear constraints (jl3p on z leading to a linearly equality 
constrained quadratic minimization problem, which can 
directly be solved by a matrix equation as in the uncon- 
strained case dSl) . The strategy is illustrated in figure [U 
which shows schematically the elliptical isolines of / for 
fixed z as a function of (zr,z.i). The values of the pa- 
rameters z determine the shape and the position of the 
parabola as well as the slope of the plane f — 0. In 
figure [U Zabs denotes the absolute minimum of /, ob- 
tained from equation ([3]). The plane / = (equation 
(ITil) ) divides the 2np-dimensional (z^, Zi)-space into the 
two half-spaces / < and / > 0. The point Zcon is the 
constrained minimum of / in the half-space / > 0, which 
lies on its boundary, i.e. on the plane / = as explained 
above. 

As long as / > /min, Zabs determines the evolution of 
the parameters. However when / = /min is reached, then 
Zcon is taken for the further integration of the trajectories 
z{t) until Zabs, driven by the constrained evolution of 
the parameters, eventually crosses the plane / = from 
/ < to / > . At this point, Zabs and Zcon coincide 
and Zabs is taken again for further integration, since / > 
leads to an increase of f{t) with time, according to the 
constraint. 

For the extension to multiple, say m, active constraints 
the real scalar valued function /(z^, z^) is simply replaced 
by the real vector valued function f(zr,Zi) = f(z) ~ 

(/!,...,/„) GR™. 




FIG. 1: The ellipses schematically represent isolines of / in 
equation ([1} for fixed parameters z. The domain of allowed z 
for the minimum of I is the full space, when / > /min and is 
reduced to the half space / > 0, when / = /min is reached. 



Now that the nonholonomic nonlinear inequality con- 
straints on z are reduced to the holonomic linear 
equality constraints (fH)) on z by the constrained TDVP, 
we can determine the constrained minimum Zcon by a 
standard method like Lagrangian multipliers. Alterna- 
tively, the constrained minimum can also be obtained 
by elimination of the dependent variational parameters. 
We prefer the method of Lagrange multipliers due to its 
generality. The method of Lagrange multipliers yields a 
compact form of the equations of motion for arbitrary 
constraints and the conditions for switching off the con- 
straints are obtained with only little additional numerical 
effort as will be shown below. Both methods however, re- 
quire a minimization problem with equality constraints. 
When inequality constraints are applied, the elimination 
of dependent variational parameters is not possible. 

We construct the function 



L 



XMz 



(15) 



with the Lagrangian multipliers A G and the real 
valued m x 2nr, matrix M = SI. The minimum of / 
under the constraint is found by dL/dui = where 




We obtain a set of linear equations 



K 




M 







with 



K 




(16) 



(17) 



where the matrix K and the vector h are the complex 
quantities of equation ([3]). If no constraint is active, i.e. 
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TO = 0, then equation (jl6p obviously reduces to the real 
formulation of equation ([3]). We use a real formulation, 
i.e. complex quantities are split into their real and imag- 
inary parts, because real constraints like / > /min natu- 
rally lead to real Lagrangian multipliers. 

The constraint (fT4|) is switched off again when Zabs 
crosses the plane / = from / < to / > 0. Finding 
this event can be accomplished in two ways. The triv- 
ial but computationally expensive way is to calculate not 
only Zcon from (jl6p . which is needed for integration, but 
additionally Zabs (from equation ([3])) after every time step 
of integration and to check when /Izaba changes its sign. 
This inefficient procedure would require the solution of a 
complex Tip X rip matrix equation for Zabs and addition- 
ally the solution of the real {2np + m) x {2np-\-m) matrix 
equation for Zcon- However it is much more efficient to 
check when A changes its sign for the special case m — 1. 
If more than one constraint is active, m > 1, it is recom- 
mended to solve the matrix equation p6p by decomposi- 
tion into two blocks, as indicated by the horizontal line 
in equation (|I6p . namely into 

Kz + M^X = h (18) 

obtained by the upper part of equation and the 

lower part 

Mi = 0, (19) 

which represents the active constraints. The solution for 
the unknowns z, A is obtained by first solving equation 
(UHl) for z 

i = K-^h~ K-^M'^X (20) 

and inserting it in equation (|19p in order to eliminate z. 
The result is a small m x m matrix equation for deter- 
mining A 

MR-^M^^ X = MK-^h e R'". (21) 

mxm 

The conditions for switching off any of the active con- 
straints are now contained in the right hand side of equa- 
tion ([2T|) . since 

f = — ^abs = M^abs = MK-^h (22) 
dz, 

due to the definitions. The ith active constraint {I < i < 
to) is to be switched off when the ith component of f Iz^ba 
changes its sign from minus to plus. 

When we insert the Lagrange multipliers calculated 
from ((2T|) in (|20|) we obtain Zcon, needed for propaga- 
tion. Numerically, the calculation of K~^\i and K~^M^ 
in (|2ip requires only one factorization of the large matrix 
K. After multiplying with M from the left the small set 
of linear equations (|2ip for determining A is obtained. 
Compared to the factorization of K the solution of the 
mxm matrix equation (|2ip for the Lagrange multipli- 
ers is negligible, since the number of parameters n will 



in general exceed the number of constraints m by far, 
e.g. in our numerical calculation there is 2np = 240 and 
the number to of simultaneously active constraints is not 
larger than three. 

B. Inequality constrained TDVP applied to GWP 

When GWP are used as trial function, it is conve- 
nient to formulate a set of linear equations for the co- 
efficients V = Vr -I- ivi first and then to obtain z from 
PH)) in a second step, just as was done in section [Til For 
these coefficients and , summarized by the notation 
(vr, Vi) = V, a similar set of linear equations is obtained. 
Equations (|10p which describe the connection between 
the time derivatives of the parameters and the coeffi- 
cients, are written in real formulation, where all complex 
quantities are split into their real and imaginary parts. 
We obtain 

A';. = -^Vl-2{{A>;y-{A'^r), 

Ak _ _l-\/k _ Ak Ak _ n Ak Ak 

^fe ^ A'^vjj + A'^Vj^q'' -t- p*", 

_^^pkfyk^k_ ^ pk/^kyk^^k _ ^j.^fe ^ l(pfc)2^ 

(23) 

with A*-- = i(4*^)-i. 

Using the notation z = (Al , , , , 7^ , , . . . , 
A'^ ,Af the complete set of equation 

(1^ for all fc = 1, . . . , iV, which are Hnear in ivQ,vf[,V2), 
may be written in short form z — U'v + d. The matrix U 
is block-diagonal with N blocks. Each block consists of 
those coefficients in equation ((23)) linear in (wq ; '^1 > ^2^)- 
The constant terms are absorbed in the vector d. The 
linear equality constraint ([14]) for a GWP trial function 
reads 

where the notation 

§^/r = E 9(if);7^^'^'^' ^^^^ 

is used. Expressing the time derivatives in equation ([^ 
by the coefficients and using ([^5)) , to arbitrary con- 
straints (f = (/i, /^) e R") imply 

f =^;7v+|Id = i7v + d = 0, (26) 
oz oz 

and hence a set of linear equations for (vr,Vi) and the 
Lagrange multipliers A e R™ is obtained 
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(27) 



with 



(28) 



Here, K — Kr + iKi and the vector r = r.^ + ivi are the 
matrix and the right hand side of equation respec- 
tively. 

We now have aU equations needed for propagation of 
coupled GWP subject to arbitrary constraints HI]). In- 
stead of ^ we solve for (v^, v^) (when no constraints 
are active both sets of equations are equivalent) after 
each time step. These coefficients are inserted in (or 
equivalently in to obtain the time derivatives of the 
Gaussian parameters, which are needed by the integra- 
tion routine to integrate the next time step. 

In order to find convenient constraints it is necessary 
to investigate the reasons for the numerical matrix singu- 
larity. The generic reasons for an ill-conditioned matrix 
K are twofold. One cause is a strong overlap of neigh- 
boring GWP, the other cause is widely spread norms of 
the GWP. A restriction on the norm of the GWP 

gmin < \ \g''\\ < g max; 

fc=l,...,iV (29) 

turns out to be sufficient to regularize the equations of 
motion. It is however more simple and numerically effi- 
cient to impose the restrictions 



fn 



7min < /''(z) = Ini7'' < 7max = fn 



(30) 



with k = 1,...,N on the amplitude of the GWP. Both 
restrictions ((29|) and ((30)) are equivalent for frozen GWP 
and they are similar even for thawed GWP (at least 
for bounded systems where the width of the GWP is 
bounded by the potential). For the active constraints 
ili = 7min or 7f = 7max) equation ([14]) using trans- 
lates into 



'•v^,-iq'=F2^q'=-ftrAj: = 0. (31) 



i^ = -vm - q' 



Therefore, in the notation of equations ([261) and (pS]) the 
entries of U are mostly zero except for the terms of equa- 
tion ([3T|) and d — tr^^. This especially simple case of 
constraints, where Gaussian parameters are bounded di- 
rectly, leads to simply temporary freezing these param- 
eters 7f when 7^ = 7,nin (7^ = 7max) is reached. As 
mentioned above, the equations of motion can instead of 
using Lagrange multipliers be alternatively obtained by 
elimination of the dependent parameters. The frozen 7*^ 
must be simply ignored in the variation. However addi- 
tional calculations are then necessary to find the criteria 
for switching off the constraints. 



Should in some cases the restriction on the amplitudes 
([50]) not be adequate, an upper bound on the maximum 
of the allowed overlap of neighboring GWP or a lower 
bound on the least eigenvalue of the matrix may be ap- 
plied. 



IV. NUMERICAL RESULTS 

Numerical tests using coupled GWP were often per- 
formed in one dimension, e.g. on the Morse potential 
0i 0, ■ Here we use a two dimensional non-integrable 
potential for testing our method. The Hamiltonian of the 
system represents the diamagnetic Kepler problem in a 
2D rotating {x,z) frame (for review, see e.g. 0, HI]). 
The magnetic field axis is directed along the z-axis. The 
potential in regularized semiparabolic coordinates reads 



with 



2 2 
r ~ X 



2 2 

z , n = r 



V = r — z. 



(32) 



(33) 



The parameters are set to a = 1/2 and /3 = 1/5 in our 
calculations. 

The method of free GWP propagation is compared to 
the method of constrained GWP propagation. The value 
of the lower bound in equation (j30p is 7min = —6.5, an 
upper bound was not needed. The comparison is pre- 
sented in figure O The trial wave function consists of 
eight GWP with the same initial values for both calcula- 
tions. Solid lines represent results of the free propagation, 
dashed lines represent the results of constrained propa- 
gation. In figure [^l^a) that normalization parameter 7^ (i) 
is selected and drawn that first reaches 7min = —6.5 at 
t K, 6.9 tc; where td is the classical period of small har- 
monic oscillations around the minimum of the potential. 
This choice allows for a direct comparison, because the 
trajectories of both calculations are equal before 7min is 
reached for the first time by any of the 7*"', and they differ 
afterwards. Normalization parameters of the other seven 
GWP are not plotted but show similar qualitative be- 
havior. In terms of figure [H the trajectories 7^ (i) in the 
range 6.9 td < 7.15 td are obtained using i^hs for the 
integration of the solid line and using Zcon for the integra- 
tion of the dashed line. Obviously the trajectory repre- 
sented by the dashed line sticks to the value 7„iin ~ —6.5 
till t « 7.15 tci where Zabs crosses the plane 7^ = 0. This 
scenario repeats several times as can be seen in the figure. 
Figure [2jb) compares the step sizes used by the variable 
step Adams routine to integrate the trajectories. The 
integration of the unconstrained equations of motion be- 
comes extremely slow around t « 7. ltd, and later on 
again for several times where the step sizes become tiny. 
Obviously there is a strong correlation between very low 
values of 7*^ in panel (a) and extremely small step sizes 
in panel (b) for unconstrained propagation. In regions 



7 





/ \ 


-4 




-6 






-8 


Ymin=-6.5 i 






(a) 1 





1e-04 



1e-06 



1e-08 




.5 9 



O 

CD 
DC 



O 



0.6 
0.4 
0.2 


-0.2 

-0.4 

0.14 
0.12 
0.1 
0.08 
0.06 
0.04 
0.02 





(b) 



4 6 

t/tcl 



8 10 



FIG. 2: Comparison between free and constrained GWP prop- 
agation of the same initial wave function consisting of a su- 
perposition of eight GWP in the 2D diamagnetic Kepler prob- 
lem: (a) Normalization parameter that first reaches 7min, (b) 
comparison of the step sizes At used by the integrator (vari- 
able step Adams method). Solid lines denote free propaga- 
tion, dashed lines denote propagation with the constraints 
7*^ > 7min = —6.5, k — 1,...,N. In both calculations the 
same error tolerances were used. 



where the free propagation is very slow, the step sizes for 
the constrained propagation are about two to four orders 
of magnitude larger, resulting in a much faster integra- 
tion. 

The magnitude of / at its minimum is a measure of 
the accuracy of the variational approximation (Tol . |20| . 
Therefore a comparison of the minima /Iz^ba ^^'^ ^kcon 
allows for an estimate of the loss of accuracy introduced 
by the constraints. A comparison of the minima shows 
that /|zco„ is slightly increased at f « 6.9 <ci with respect 
to /|zaba but at later times, one approximation is about as 
good as the other in the average, although a poorer ap- 
proximation of the constrained wave function to the exact 
one would be expected. However it has been shown that 
the approximate wave function determined by TDVP is 
not always the 'best' possible approximation of the trial 
function to the exact wave function [20]. There might 
be regions on the manifold of the trial function that are 
closer to the exact wave function than the function deter- 
mined variationally, especially when the manifold has a 
large curvature and long time intervals are considered. 
This fact, together with the insensitivity of the wave 
function to small variations of the parameters in some 



FIG. 3: (a) Real part of the auto-correlation function. The 
initial wave function for the 2D diamagnetic Kepler problem 
is a superposition of 20 GWP. Variational propagation with 
the constraints 7* > —6.5 (dashed line) is compared with nu- 
merically exact calculations (solid line) and with frozen GWP 
propagation (dotted line). The results of the constrained 
calculation and the numerically exact calculation practically 
coincide, and nearly no deviation is visible, (b) Deviations 
of the auto-correlation functions calculated by constrained 
GWP (dashed line) and by frozen GWP (dotted line) from 
the exact one. 



directions in case of a singular matrix, may explain the 
behavior of only temporary slight loss of accuracy intro- 
duced by the constraints. The insensitivity of the trial 
wave function to the constraints can also be deduced from 
the auto-correlation functions C{t) = {x{t = 0)\x{t)) ob- 
tained by both methods since they almost coincide and 
no deviation from each other could be seen in any figure. 

To demonstrate the accuracy of the constrained GWP 
method a superposition of 20 GWP having all the same 
width and zero momenta, equally distributed on an 
equidistant grid was used as the initial wave packet. 
This initial wave packet was propagated by three dif- 
ferent methods. The real parts of the resulting auto- 
correlation functions are plotted in figure[3lja) . The imag- 
inary parts, not shown in a the figure, exhibit similar be- 
havior. For reference the numerically exact propagation 
was performed by the split operator method [21| (solid 
line). The result of our constrained (7*^ > 6.5) GWP 
propagation (dashed line) is mostly very accurate and 
nearly no deviation from the exact solution is visible for 
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many classical periods. By contrast, the result obtained 
from a frozen Gaussian propagation (dotted line) turns 
out to be much more inaccurate. This becomes partic- 
ularly apparent in figure E^b) , where the deviation be- 
tween the exact time signal and the time signals obtained 
from constrained (dashed) and from frozen width (dot- 
ted) propagation is plotted. For short times both meth- 
ods are very accurate and nearly no deviation between 
the time signals is visible. With increasing time, however, 
the accuracy of the frozen width calculation is lost much 
faster than that of the constrained propagation. This 
is not completely unexpected since the constrained trial 
function still has more free variational parameters than 
the frozen GWP method and therefore the constrained 
calculation is slower. Note that an unconstrained propa- 
gation of these 20 GWP with variable widths according 
to the TDVP would not be possible. With the propa- 
gated wave packet at hand it is straightforward to obtain 
e.g. the eigenvalues of the Hamiltonian by Fourier trans- 
form or harmonic inversion (for a review see e.g. [23 |) of 
the auto-correlation function or to extract eigenfunctions 
of the system (e.g. [23j). 



has been proposed. The method is based on applying 
nonholonomic inequality constraints on the motion of 
the GWP. The constraints must be chosen to prevent 
the matrix from becoming singular. From the inequality 
constrained TDVP a simple matrix equation for the time 
derivatives of the parameters is obtained just as in uncon- 
strained TDVP. The method is in fact applicable for ar- 
bitrary trial functions and inequality constraints. For the 
GWP trial functions we found it sufhcient in most cases 
to apply simple bounds on the normalization parame- 
ters to regularize the matrix and to obtain well-behaved 
equations of motion, rendering the integration orders of 
magnitude faster. The loss of accuracy of the method 
caused by the constraints is found to be negligible for 
sufhciently many GWP. The method allows for the prop- 
agation of a large number of coupled GWP, as compared 
to the unconstrained GWP propagation, and guarantees 
accurate results within reasonable time. Our method for 
time propagation of constrained coupled Gaussian wave 
packets presented in this paper will be very powerful in a 
large variety of future applications to overcome the prob- 
lems with ill-conditioned and stuck differential equations. 



V. SUMMARY 



A novel method to overcome the matrix singularity 
problem in the variational Gaussian wave packet method 
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